% This script calculates PNJunction of parabolic bands with MoS2
% parameters. Various Dopant Concentrations are considered.

% Physical Constant Setup (In SI Unit)
a0 = 0.5291772109217;  % Bohr Radius in Angstrom 
eps = 1;               % Relative permittivity
eps_0 = 8.8541878176E-12;   % Vacuum permittivity 
e = 1.6E-19; % Electronic Charge
hbar = 1;
Vh = 27.2113850560;
kB = 8.6173324E-5/27.21138386;
vf = 1E6/(2.187691263373E6);

% Parabolic DOS Setup
mv = 0.58;
mc = 0.48;
Eg = 1.85/Vh;
Ev = -Eg/2;
Ec = Eg/2;
deg = 6;

% Charge Setup
f1 = -1E-6;
f2 = 1E-6;

% Spatial Grid (In Hartree Atomic Unit)
L = 500000;
Nx = 100000;
ncx = 1+Nx/2;
x = linspace(-L/2,L/2,Nx);


%% First Step
[Va, Ra, Err] = PNParab(f1,f2,mv,mc,Ev,Ec,deg,300*kB,2E-3,Nx,L);
filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
dlmwrite([filename, 'V.txt'],[x',Va']);
dlmwrite([filename, 'R.txt'],[x',Ra']);
dlmwrite([filename, 'Err.txt'], Err);


%% Second Step
filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
V0 = dlmread([filename, 'V.txt']);
Vx = V0(:,2)';
[Va, Ra, Err] = PNParab(f1,f2,mv,mc,Ev,Ec,deg,300*kB,2E-3,Nx,L,Vx);
    
filename = sprintf('MoS2(f1=%1.1E,f2=%1.1E,L=%i,N=%i)',f1,f2,L,Nx);
dlmwrite([filename, 'V.txt'],[x',Va']);
dlmwrite([filename, 'R.txt'],[x',Ra']);
dlmwrite([filename, 'Err.txt'], Err);

%% Plotting data
L = [1000,2000,3000,4000,5000];
N = [2001,4001,6001,8001,10001];
linespec = {'-b', '-c', '-m','-y', '-g'};
figure; hold on;
%ax1 = subplot(2,1,1); xlim([0,200]); hold on;
ax2 = subplot(1,1,1); xlim([1,500]); hold on;
for i = 1:5
    filename = sprintf('MoS2(f=%1.1E,L=%i,N=%i)',1E-4,L(i),N(i));
    V = dlmread([filename, 'V.txt']);
    R = dlmread([filename, 'R.txt']);
    %plot(ax1,V(:,1),V(:,2), linespec{i}); 
    
    plot(ax2,R(:,1),1./R(:,2), linespec{i});
    
end
x = linspace(-1500,1500,6001);
R_metal = 2*(0.0318)/(2*pi^2)./x;
plot(ax2,x,1./R_metal,'-k', 'linewidth',2);
    
    
    
    
    